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Falling liquid films become unstable due to inertial effects when the fluid layer is sufficiently thick 
or the slope sufficiently steep. This free surface flow of a single fluid layer has industrial applications 
including coating and heat transfer, which benefit from smooth and wavy interfaces, respectively. 
Here we discuss how the dynamics of the system are altered by feedback controls based on observa¬ 
tions of the interface height, and supplied to the system via the perpendicular injection and suction 
of fluid through the wall. In this study, we model the system using both Benney and weighted- 
residual models that account for the fluid injection through the wall. We find that feedback using 
injection and suction is a remarkably effective control mechanism: the controls can be used to drive 
the system towards arbitrary steady states and travelling waves, and the qualitative effects are in¬ 
dependent of the details of the flow modelling. Furthermore, we show that the system can still be 
successfully controlled when the feedback is applied via a set of localised actuators and only a small 
number of system observations are available, and that this is possible using both static (where the 
controls are based on only the most recent set of observations) and dynamic (where the controls are 
based on an approximation of the system which evolves over time) control schemes. This study thus 
provides a solid theoretical foundation for future experimental realisations of the control of falling 
liquid films. 


I. INTRODUCTION 

The flow of a thin liquid film down an inclined planar wall is a classical problem in fluid mechanics. 
The flow becomes unstable when the Reynolds number, defined on the undisturbed interface flow speed, is 
above a critical value which depends on the inclination angle; the flow is stable when the layer is sufficiently 
thin. After the onset of instability, the system initially exhibits two-dimensional ( 2 -D) waves that propagate 
down the slope, followed by more complicated behaviour that can eventually lead to three-dimensional ( 3 -D) 
spatiotemporal chaos. The development of thin film models for this system, and the behaviour exhibited 
therein, has recently been reviewed by Craster and Matar^i^, Kalliadasis et al^. 

In addition to acting as a paradigm for understanding transitions between different types of dynamical 
behaviour, the flow of thin films has a broad range of industrial applications. We note particularly coating 
flowPI, where a uniform coating of a flat or shaped substrate is desired, and heat and mass transfer, which 
is typically enhanced by mixing associated with interfacial wave^. The film dynamics and stability can be 
influenced by feedback control, which may be able to delay the onset of instability, and in an ideal situation 
we would like to be able to drive the system into the full range of regimes. 

Much of the thin-film literature focuses on the additional instabilities and flow modes that can occur in 
flows with heating and cooling, or on flows over non-uniform topographji^®^, both ha ve dire ct applications 
in heat exchangers. Thermal effects and topography are often combined with each otheif^^^^ or with electric 
fieldpHni. 

Other physical mechanisms that have been investigated within the context of thin-film flow 
down inclin ed planes include chemical coatings or microstructure to induce effective slipP^, surfactantJ^, 
porou^2012ij gj. deformabl^^ walls, explicit injection/suction through the planar walP^and magnetic field^^. 
For flat homogeneous walls, a steady uniform flow solution exists known as the Nusselt solution, which can of 
course be unstable. In contrast, imposed heterogeneity can be used to create patterned states with different 
stability properties to the corresponding homogeneous system. However, systems in which the controls are 
able to actively respond to the instantaneous system state, rather than being pre-determined (open-loop 
control), have a much stronger effect on stability. In this paper, we consider such closed-loop (or feedback) 
control, which will be imposed to the system by suction and injection through the wall. 

Thompson, Tseluiko, and Papageorgioul^ used long wave models to study the effect of imposed, steady, 
spatially-periodic suction/injection on thin-film flow down an inclined plane. They found that the imposed 
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suction always leads to non-uniform states, enables a non-trivial bifurcation structure and complicated time- 
dependent behaviour, and significantly alters the trajectories of particles in the flow, but has a relatively 
small effect on flow stability. Fluid injection through slots has also been considered theoretically for its effect 
on spreading film^^ suction leads to ridges on the free surface, and injection to indentations, but there is 
no steady state as the total mass is not conserved. Injection has some similarities to flow over a porous 
wall, which tends to wick fluid into narrow pores; this flow is particularly relevant to the printing of ink 
onto paper, for which substrate porosity affects the lifetime and spreading of drops. Davis and Hockingl^ 
considered flow of thin drops and films with wetting fronts along a porous substrate that is wetted by the 
fluid, and is initially dry, and found that for both films and drops, the fluid is eventually drawn entirely 
into the substrate. Thiele, Goyeau, and Velarde^ used a Benney equation to study flow over a heated, 
fluid-filled, inclined porous substrate, bounded below by a solid wall so that the total mass of the liquid film 
is conserved. They found that the addition of a porous substrate typically has a small destabilising effect on 
the uniform film state, and in the nonlinear regime the film develops drops and ridges which slide down the 
plane. Film flows with point and slot-shaped sources have also been studied as a model for lava spreadin^^, 
though the fluid does not form a continuous film, instead containing several wetting fronts. 

Feedback control requires observations of at least some components of the ^tem state, and we will build 
our control strategies around observations of the film height. Liu and Gollub^ investigated experimentally 
the dynamics of thin films within the context of the onset of chaos; they used a fluorescence imaging process 
to measure the two-dimensional film thickness in real time, and also used laser beam deflection to obtain 
local measurements of the interface slope. Vlachogiannis and Bontozoglou^ examined the flow of thin 
films over a wavy wall, and used interferometry calibrated against needle-point measurements to obtain 
the interface height. Heining, Poliak, and Sellier^^ showed that the free surface shape and topography 
profile can be obtained from measurements of the surface velocr^ and implemented this both in Navier- 
Stokes simulations and experiments. Schorner, Reck, and Aksel^ used experiments with visualization by 
laser reflection to study the effect of differently shaped topographical configurations with the same basic 
amplitude and wavelength on the flow down an inclined plane; they were able to infer the streamwise 
growth rate of small-amplitude perturbations by comparing the magnitude of interfacial fluctuations at two 
streamwise locations. 

Thin film flow is often studied using reduced-dimensional models, which differ most fundamentally in the 
manner in which inertial effects are incorporated. Here we use two different first-order long-wave models: the 
Benney equation and the weighted-residual (WR) equation, which were extended by Thompson, Tseluiko, 
and Papageorgioul^ to include the effect of suction and injection. These two long-wave models are identical 
at zero Reynolds number, and both agree with the Navier-Stokes system regarding the critical Reynolds 
number for the onset of instability. However, the structures of these models differ significantly, in particular 
the number of degrees of freedom. The Benney equation is a single evolution equation for the interface height 
h(x,t), while the weighted-residual model comprises coupled equations for and the independently- 

evolving down-slope flux q(x,t), and of course the Navier-Stokes equations at finite Reynolds number allow 
evolution of h(x, t) together with evolution of the flow velocity at every point within the fluid. The robustness 
of control strategies to changes in the model is one of the major themes of this paper; we seek to understand 
what features of the system state need to be measured to deliver effective control, and whether the control 
system can be designed without needing detailed knowledge of the system state and underlying dynamics. 

Feedback control systems consist of a set of control actuators and response function^^ for linear controls 
the response is a linear function of the deviation of the observed state from the desired state. An appropriate 
linear function is constructed based on hypotheses regarding the dynamics of the uncontrolled system and 
its response to the control actuators. In the simplest scenario, the controls are distributed along the domain. 
In more realistic scenarios, the control actuation is only possible at a finite number of points in the domain, 
and observations of the current state are not available everywhere. Feedback control theory is useful in both 
cases, and standard tools for tackling these problems are discussed by Zabczyk^^. 

Theoretical applications of feedback control for thin films include thermal perturbations in liquids spread¬ 
ing over a solid substrate to suppress the contact line instabilitjEH and point actuated suction/i^ction to 
suppress waves in weakly nonlinear models of thin film^^ or to enhance such wavy behavioui®. When 
suppressing waves, again in weakly nonlinear systems, Armaou and Ghristofides^ prove that stabilisation 
can be a chieved by using a finite number of observations, either by static or dynamic output feedback control 
(see Sec. IVD and appendix]^, while Armaou and Ghristofidesl^ use nonlinear feedback controls. Efforts 


have also been made to optimise the placement of actuators and sensors to suppresJ^ or enhancd^ waves 
on a thin film. Feedback control strategies have been implemented for the two-dimensional Navier-Stokes 
equations in the context of data assimilatioiP^, in which controls applied towards known observations are 
used to overcome incomplete knowledge of the initial state in the forecasting of hurricanes and typhoons. 
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Both of the long wave models studied here reduce to the forced Kuramoto-Sivashinsky (KS) equation 
under a we akly nonlinear analysis, and form part of a rational hierarchy of models which lead to the 
KS equatiorPSEni. xhe KS equation retains many essential features of the Benney equation, in particular 
nonlinearity, energy production and dissipation. However, it is a much simpler system, which moreover 
has the property of global existence of solutionPil, a global attractoi!^, and sharp bounds on its solutions 
and derivatives^. These global properties make the KS equation amenable to analysis, and various control 
schemes have been developed and applied for the KS equation, e.g. Gomes, Papageorgiou, and Pavliotisl^. 
With the inclusion of suitable linear feedback controls towards the flat solution, the bounds on the solutions 
make it possible to prove that the L^-norm is a Lyapunov function for the dynamical system corresponding 
to the KS equation, and hence the controlled system is nonlinearly stable. Furthermore the existence of 
bounds on the solutions can be used to prove the existence of optimal controls. When stabilising nontrivial 
steady states or travelling waves, the boundedness of these nontrivial solutions play a crucial role in the 
proof of existence of a Lyapunov function and also influence the choice of the controls. 

The difficulty in applying the KS controls to thin-film systems is twofold: the nonlinearities are more 
complicat ed, and so there is no global existence theory; in fact the Benney model can lead to unbounded 
behaviour^®!, and secondly, the structure of the weighted-residual and Navier-Stokes models, when applied 
at finite wavelength, is significantly different to that of the KS equation. It is reasonable to suspect that a 
control strategy carefully optimised for one model may be ineffective in another, and so we focus here on the 
use of relatively simple control schemes, and investigate their robustness to variations in the model details. 


In this paper, we consider the effect of feedback control, applied via a suction boundary condition and 
based on observations of the interface height, on the dynamics and stability of a thin layer of fluid flowing 
down an inclined plane. We use two different long-wave models, described in Sec. to model the system, 
and also compare certain stability results to those for the Navier-Stokes equations. In Sec. |III[ we show 
that a proportional control scheme, in which fluid is injected at each streamwise location in proportion 
to the observed deviation of the interface height at that location from uniform, has a stabilising effect on 
nearly-uniform flow in all three models. We compute the critical control magnitude required to stabilise 
the uniform state to perturbations of all wavelengths, finding that the critical Reynolds number for the 
onset of propagating waves can be increased significantly by using the proportional control scheme. For the 
control system to be physically realisable, we must relax the requirement that the system state is known at 
every position, and that the suction applied can take an arbitrary shape. Therefore, in Sec. IV we discuss 
control strategies for the case where feedback is applied through a discrete set of localised actuators and the 
system is similarly observed at a small number of locations in the domain. We discuss both control schemes 
based on static observations (where the control amplitudes are calculated from only the most recent set of 
observations) and dynamic observations (where the controls are calculated using an approximation to the 
system state which is built up over time). In Sec. we discuss linear stability and nonlinear behaviour 
when controlling to non-uniform travelling waves and non-uniform steady states. We note that the property 
of a given non-uniform state being an exact solution of the equations is model dependent and therefore 
can never be perfectly satisfied. In order to test the robustness of the control scheme to variations in the 
model, we consider controlling to non-uniform interface shapes with only very crude prior knowledge of 
the governing equations. We find that in this case, the controls lead to equilibrium states which converge 
towards the desired system as the control amplitude is increased, and that these equilibrium states are stable 


when large-amplitude controls are applied. Our conclusions are presented in Sec. VI 


II. GOVERNING EQUATIONS 


We consider a thin layer of fluid, with mean thickness hg, flowing down a plane inclined at an angle 9 
to the horizontal. We adopt a coordinate system such that x is the down-slope coordinate, and y is the 
perpendicular distance from the wall, as shown in Fig. The upper interface of the fluid is a free surface, 
located at y = h(x,t). The lower boundary of the fluid is a rigid wall, through which fluid may be injected 
or removed. 

The two-dimensional (2-D) Navier-Stokes equations admit a solution which is uniform in the streamwise 
direction, known as the Nusselt solutioiP^, for which the surface velocity is Ug = pgsin 9/{2y), where p 
is the fluid density, g the acceleration due to gravity and p the dynamic viscosity of the fluid. We non- 
dimensionalise the problem based on the length scale hg and the velocity scale Ug, and define the Reynolds 
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FIG. 1. Sketch of flow domain showing coordinate system. We consider a fluid layer, with mean height hs, bounded 
on y = 0 by a rigid wall inclined at an angle 6 to the horizontal, and at i/ = h{x) by a free surface. Fluid is injected 
through the wall, with velocity v = F{x,t) which changes in time in response to fluctuations of the free surface. 


number R and the capillary number C based on the velocity Us'. 

n phsUs pXJs 

it - ? ^ - 5 

M 7 


( 1 ) 


where 7 is the coefficient of surface tension at the air-fluid interface. Subsequent equations are all dimen¬ 
sionless. 


A. Navier-Stokes equations 

We wish to solve the 2-D Navier-Stokes equations, with velocity u{x,y,t) = {u,v), and fluid pressure 
p{x,y,t). The momentum equation and continuity equations are 

R {ut + UUx + VUy) = -Px + 2 + Uxx + Uyy, (2) 


and 


R {Vt + UVx + VVy) = -Py - 2 cot 6 > -b Vxx + Vyy 


Ux+Vy = 0 . 

The boundary conditions at the wall are given by 


( 3 ) 

( 4 ) 


u = 0, V = F{x, t). 


( 5 ) 


Here the function F{x, t) represents the injection velocity normal to the wall, y = 0. Note that we assume 
that the injection of fluid does not affect the no-slip boundary condition on the wall. At the interface, 
y = h{x,t)i the tangential and normal components of the dynamic stress balance condition yield 


P-Pa 


F tty) hx') T 2hx {Vy Ux) - 0, 

2 

2 _|_ ^2 F Uxhx — hx {Vx F Uy)) = ~ 


1 hxx 

C(1 + /i2)3/2’ 


( 6 ) 

( 7 ) 


where Pa is the atmospheric pressure, assumed constant. The system is closed by the kinematic boundary 
condition at the free surface 


ht = V — uhx 


( 8 ) 
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Defining the down slope flux q 


qix,t) = / u{x,y,t)dy, (9) 

^0 

integrating Q, and applying the boundary conditions ([^ and (|^, yields the mass conservation equation in 
terms of q: 

ht - F{x,t) + q^ = 0. (10) 

The flow is modelled either by the 2-D Navier-Stokes equations, or by one of two reduced-dimension 
long-wave models, which are derived according to either the Bennej^ll or weighted-residual methodology 
in order to approximate the flow in a long wave limit. To achieve this, we define new variables 

X = 5x, T = 6t, V = 6w^ C = (5^(7, F = Sf (11) 


and seek a solution for small <5. The extension of these two long wave models to include the effects of 
suction/injection is discussed by Thompson, Tseluiko, and Papageorgioul^, and so the relevant governing 
equations are stated without derivation here. 


B. Benney system 


In the Benney model, the flux q is slaved to the interface height h, and up to first order in 6, including 
the cross-flow effects induced by F, the flux is given bji^ 


2h^ 

q{X,T) = ^+8 


— ( —2hx cot 9 + 


hxxx \ ^ ^ Sh^hx 


C 


-J 


V 15 


We then recast the equation in terms of the original variables to obtain 


q{x, t) = — [2 — 2hx cot 9 + 


+r(^ 


Sh^hx 2h‘^F 


= Z{h,F). 


( 12 ) 


(13) 


The coupling of (13) to (10) yields a closed system for the evolution of the interface height h{x,t). 

We note that the appearance of terms involving F in (13) is a consequence of the choice of F with respect 
to the long wave scaling. By supposing F to be an order smaller with respect to 5 in the long wave expansion 
(11), we can replace (13) with the simpler version: 


q{x, t) = — ( 2 — 2hx cot 9 + 


C 


SRh^hx 

1 ^ 


(14) 


In this limit, the only effect of F on the system dynamics is via the mass conservation equation (10). 


C. Weighted-residual system 


Alternatively, following the weighted-residual methodology developed by Ruyer-Quil and Mannevillel^, 
the flux q gains its own evolution equation, so that time derivatives of both h and q appear in the equations. 
After substituting © into the governing equations and retaining terms up to and including 0{S), the 
evolution equation for q i^^ 


26Rh‘^ dq 2h^ , 

— 3? + ’=-+'* 




R 


fISq^hx Sdhqqx 


35 


35 




In the original variables, (15) becomes 
2 


-Rh^qt +q = — 
5 3 


2 — 2hx cot 9 + 


C 


R 


/ ISq'^hx Sdhqqx hqF 

3^ ^ 


= Z{h,q,F), 


(15) 


(16) 


which when coupled to (10) yields a closed system, for h{x,t) and q{x,t). Initial conditions are required for 
both h and q. The Benney and weighted-residual models are identical when i? = 0, and can be shown to 
agree at 0(1) and 0{5) in the long-wave limiiP^. 
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D. Choice of controls 

The focus of this paper is on the application of suction as a linear control mechanism in response to 
observations of the interface height. We begin in Sec. |III| by considering the case of controlling towards the 
uniform Nusselt state, based only on observations of h. To achieve this, we set 

F{x,t) =—a[h{x,t) — 1], (17) 

where a is a real constant to be chosen; in most cases we find that the uniform state becomes increasingly 
stable for large positive a. Note that if h = 1 everywhere, then the controls have zero magnitude. 

The control scheme 0 requires perfect knowledge of the instantaneous interface shape h(x,t), and the 
ability to impose any continuous F{x^ t). In practice, we expect neither of these assumptions to hold. Instead, 
fluid is injected via a number of localised actuators, or slots, in the substrate, and interface observations 
are available at a small number of locations in the flow domain. In Sec. |IV[ we investigate control schemes 
based on point actuators and localised observers, with both static observations and dynamic observers, and 
the appropriate form of the injection profile F{x,t) will be discussed when necessary. 

In Sec. |V] we consider controlling towards either nonuniform travelling waves of permanent form or 
nonuniform steady states. Travelling waves can be written as h = i7(C), where (^ = x — Ut and U is the 
constant propagation speed. By direct analogy to 0 , we set 

F(C,t) = -a[h(C,t)-i7(C)]. (18) 

We note that if h{x, t) = H{x — Ut) for all time, then f = 0, so that the travelling wave h = H{x — Ut), is 
also a solution of the controlled equations. Nonuniform steady interface shapes H{x) are not steady states 
of the equations when F = 0, but Thompson, Tseluiko, and Papageorgiou^^ showed that imposing a steady 
suction component S{x) enables non-uniform steady states. Combining with linear control, we obtain 

F(a;, t) = —a[h{x, t) — H{x)] + S{x). (19) 

For non-uniform states, the calculation of S{x) to obtain an exact steady state, or of the travelling wave 
solution H{(), require detailed knowledge of the governing equations. For example, these states differ even 
between the Benney and weighted-residual models, let alone the Navier-Stokes equations. In Sec. |V C[ we 
consider the robustness of our control schemes when the model details are not well known; we do so by 
controlling towards a finite-amplitude non-uniform state H{x), but setting S{x) = 0, so that the target 
state is not a steady solution. As a result, the control parameter a has a role to play in setting both the 
shape of any steady states obtained, as well as their stability. 


E. Numerical calculations for linear stability of non-uniform solutions 


For a translationally invariant system, as occurs for distributed controls towards a uniform him state, 
the linear stability of the uniform him state can be calculated via a normal mode analysis, and this will be 
pursued for the Benney, weighted-residual and Navier-Stokes systems in section |III[ However, if the base 
state for the stability analysis is not uniform, or the feedback control system has localised actuators or 
observers, then the system is no longer translationally invariant, and so the eigenmodes of the system are no 
longer normal modes. In that case, we can compute the discretised eigenmodes of the system by formulation 
and numerical solution of a generalised eigenvalue problem for linear stability, as described below. 

We consider the evolution of a small perturbation h: 

h = F[ (x) + ehe^*, q = Q{x) + eqe^^, F = S{x) — ee^^ah. (20) 

We recall the Benney equation: 


ht+qx-F = 0, q = Z{h,F), 


( 21 ) 


where Z{h,F) is dehned in (13), and expand for small e. The equations at 0(1) in e must be satished by 
the base state H{x), Q(x), S{x). At 0(e), we obtain a generalised eigenvalue problem for h, q and A: 


A 




f -al -dA fh\ 
{Zh - Zpal -I )\qj 


( 22 ) 
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where dx is the derivative operator, I is the identity matrix, and the blocks Zh and Zp are linear operators, 
for example from (13) we have 


Zh = 


2 - 2 cot OHx + 


Hx 


C 


IQRH^Hx SH^RS 


2H^ 


cot 6 


15 


a -I_ a 


We can eliminate q to obtain a smaller eigenvalue problem for h alone: 

Xh = —dxZhh — [/ — dxZp]ah. 


(23) 


(24) 


For the uniform state, H = 1, Q = 2/3, and S' = 0, the blocks Zh and Zp simplify considerably; in fact, 
we can calculate the eigenvalues analytically in that case. For non-uniform base states, we calculate the 
eigenvalues by replacing the derivative operators with pseudo-spectral derivative matrices, and solving the 
eigenvalue problem numerically using standard algorithms available in Matlab. 

We can also write the flux equation of the weighted-residual system in a similar form: 


‘^Rh'^qt +q = Z{h,q,F). 

We again obtain a generalised eigenvalue problem for /i, q and A in the weighted-residual equations: 


0 iRH'^I 


—al —da 
Zh — ZpOtl Zq — I 


(25) 


(26) 


where the blocks Zh, Zq and Zp can be obtained by differentiating (16). We note that there are twice as 
many eigenmodes in the weighted residual equations as in the Benney equation. 

In addition to linear stability calculations, we also perform a number of initial value calculations as 
described by Thompson, Tseluiko, and Papageorgioul^, involving a pseudo-spectral method for spatial 
discretization, and a fully-implicit, backward finite difference time stepper. 


III. EFFECT OF PROPORTIONAL CONTROLS ON THE STABILITY OF A UNIFORM FILM 

The uniform film state h = 1, known as the Nusselt solution, is a steady solution to all three sets of 
equations (Navier-Stokes, Benney and weighted-residual) in the absence of suction. The base state is 

h = l, (? = 2/3, u = y{2-y), v = 0, p = 2(1 - y) cot 6». (27) 

In 2-D Navier-Stoke^^^*^, Bennej® and weighted-residual model^^, this solution is linearly stable to per¬ 
turbations of all wavelengths if 


5 

R < Rq = - cot 9. (28) 

As R is increased across this threshold, the first perturbations to become unstable are those with infinite 
wavelength, and in fact the long-wavelength nature of the instability was the physical motivation for the 
development of long-wave models. 

The application of linear controls F = —a{h — 1) affects the linear stability of the Nusselt solution. As 
the system is invariant under translation in x, the eigenmodes are proportional to exp{ikx), and so we write 

h = 1 + ^-f (29) 

3 


and seek a solution for e ^ 1. We aim to compute X{k)-, solutions are stable to perturbations of all 
wavelengths if the real part of A, 3?(A), is negative for all real k. In what follows we calculate A for each of 
the models including feedback control in order to establish that the constant a can be chosen to stabilise the 
uniform flow (27). In the case of Benney and weighted-residual models this can be performed analytically, 
whereas for the Navier-Stokes equations we compute the eigenvalues numerically. 
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FIG. 2. Results for the real part of the Benney eigenvalue A as a function of k, for R = 5, C = 0.05, 6 — 7r/4 and 
a = 0 (solid lines), and a = as (dashed lines) from (34l. 


A. Benney equations 


The linearised mass conservation equation (101 yields 


\h + ah + ikq = 0. 


Substituting (29) into (13) gives 


q= 2- 


2ik cot 9 


ik^ 


8ikR 2aR' 


15 


3 J 


and combining (31) with (30) yields a single eigenvalue A: 

2Rik\ 


A = —a 1 + 


r, 7 / 5cot0 


8C 


(30) 


(31) 


(32) 


Throughout this paper, we will suppose that a is real and independent of k, and taking a > 0 in (32) is 
seen to have a stabilising effect on the Benney system, li R < Rq, the Nusselt solution is linearly stable for 
all real k in the absence of controls, and becomes more so as a increases. However, if R > Rq, there is a 
finite k with maximum growth rate, and it is easy to show that 


max3fi(A) = —a + 

k 


16C{R-Rq)^ 
^ 75 


(33) 


Hence we can stabilise the uniform film state against perturbations of all wavelengths by choosing a > as, 
where 


OLB = 


16(7(7? - i?o)^ 
75 


(34) 


The dispersion relation (32) is plotted with and without controls in Fig. for parameters at which the 
uncontrolled solution is unstable. In the absence of controls, the real part of the A is positive for small fc, 
with a finite cutoff wavenumber fc, above which the real part rapidly becomes increasingly negative. Setting 
a = as shifts the real part of the entire spectrum by —as- This means that perturbations of very small 
wavenumber decay with a finite growth rate of approximately —as, rather than having a small positive 
growth rate in the absence of controls. The maximum growth rate occurs at the same k as in the absence 
of controls, and for a = as, this maximum growth rate is exactly zero. We can also compare the imaginary 
part of A (not shown); we Hnd that setting a = as slightly increases the magnitude of the imaginary part, 
and hence the downstream propagation speed of small perturbations is slightly increased. 
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FIG. 3. Results for the real part of both eigenvalues A for the weighted-residual model as a function of k, for R = 5, 
C = 0.05, 6 = 7r/4 and a = Q (solid lines), and a = as (dashed lines) from (34 1 . 


B. Weighted residual equations 


The linearised version of the weighted residual equations (16) yields 

8ikR 


2Ai?/ 2ikcot6 
— q + q= 2--- 


ik^ 

3C 


35 


2Ra\ ~ 68ikR „ 
h - q. 




105 


We combine (35) with (30) to obtain a quadratic equation for A: 

68ikR 2aR\ f ^ \8ikR\ 


2RX^ 


+ A 1 + 


105 


-) 


+ a 1 + 


35 


-) 


2ik - 


8k^RH Sk'^R 


15 


35 


= 0 , 


where 


5 5P 5/c^ 

R^_-cote+ — =Ro + —. 


(35) 


(36) 


(37) 


The characteristic equation (36) has complex coefficients, and so its two roots for A are not complex conju¬ 
gates. We calculate the two roots for A numerically to determine the effect of imposing controls; Fig. |^shows 
A as a function of fc, with and without controls. The eigenvalues of the weighted-residual equation display 
relatively little variation with respect to k in comparison to the Benney results, but the two systems share 
the same cutoff wavenumber in the absence of feedback controls. With the addition of feedback controls, 
we hnd that positive a decreases the real part of A for both eigenvalues of the weighted-residudal system, 
with the exception of the most stable eigenmode at fc = 0, which is independent of a (see below). Choosing 
the critical a = as for the Benney equation, given by (34), is more than sufficient to stabilise the uniform 


state against perturbations of all wavenumbers in the weighted-residual equations. 
In the long-wave limit fc ^ 1, (36) becomes 


(A + a) ( 1 + ^ 1 = 0 


(38) 


which has roots at A = —a and A = —5/(277). Choosing non-zero a affects the stability of the first root, 
and means that we must choose a > 0 to obtain a stable solution. The second root is unaffected by a, and 
as a consequence, the maximum real part of A across all k is always greater than —5/(277), regardless of the 
value of a. 

Although the effect of a on A is more complicated than that for the Benney equations, we can still 
calculate the critical control amplitude a needed to ensure that 3fi(A) < 0 for all k. Perturbations with very 
large wavenumber are always stabilised by surface tension, so have negative real part. If the uniform state 
is unstable to perturbations for some k, then there is at least one cutoff value of k for which 3f7(A) = 0. 
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We therefore investigate the conditions for which there is a purely imaginary root, writing for convenience 
A = —2ikil. We solve the imaginary part of (36) to obtain 


1 + 


O = 


9aR 


1 + 


2aR ' 


(39) 


Since is independent of k, we can rewrite the real part of (36) as a quadratic equation in k^\ 

8Rn^ 


,2 I 


136i?fl 8R 8Ro 
105 lh“ 


+ a = 0. 


(40) 


The roots of this equation correspond to wavenumbers where 3fi(A(fc)) = 0. When a is insufficient to stabilise 
perturbations of all wavelengths, there are two roots for fc^, and one root at the critical value of a. The 
uniform state is stable to perturbations of all wavelengths if there are no real roots for fc^, i.e. when (40) has 


negative determinant. This condition can be rewritten using the definition of to obtain that the uniform 
state is stable if 


R 


71aR ia^R? 


1-f 

245 

175 

1-k 

4ai? 

4a2i?2 


H-r— 


— Ro 


< 


75a 

l6C' 


(41) 


The term in square brackets is monotonically decreasing in aR for aR > 0. When a is small, we find 


R ~ Rq + 


75a 

Tec’ 


which is exactly the Benney result. At large a. 



Ro + 



(42) 


(43) 


and so the maximum R for which the uniform solution is stable at fixed a is increased by a factor of 
approximately 10 according to the weighted-residual model. The stability boundaries for the Benney and 
weighted residual results are given by (34) and (41), and are plotted in the ^/a-R plane in Fig. together 
with the corresponding Navier-Stokes results as discussed below; it appears that the stable region of the 
Benney equation is always a subset of the stable region according to the weighted-residual equation, so the 
critical a predicted by (34) is indeed a conservative estimate of the necessary a required to stabilise the 


uniform film to perturbations of all wavelengths. 


C. Navier-Stokes equations 


We can compute the linear stability of the Nusselt state in the two-dimensional Navier-Stokes equations, 
subject to distributed feedback controls, by a normal mode analysis. This analysis is well known in the 
absence of suctioiP^. The addition of suction controls changes only one boundary condition in the resulting 
Orr-Sommerfeld system, and so only a brief description of the equations is presented here. 

We perturb about the uniform state, writing 


h = 1 
u = u{y) 

V = 9 

P = P{y) 


eH exp(ifca; -I- At) 
eU{y) e'xjp{ikx + At) 
eV (y) exp{ikx + At) 
eP{y) exp{ikx + At), 


(44) 


where u{y) = y(2 — y) and p{y) correspond to the uniform film solution described in (27), and then linearise 
with respect to e. The perturbation velocity components tJ{y) and V{y) can be expressed in terms of a 
streamfunction '4>{y), so that 


U{y) = -^'{y), V{y) =iki>{y), 


(45) 
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FIG. 4. The boundaries for stability to perturbations of all wavelengths, for 6 = tt jA, C = 0.05. The stable region 
emanates from the axis. 


which immediately satisfies the mass conservation equation The two components of the momentum 
equation ([^ and ([^ can then be combined to yield the Orr-Sommerfeld equation, which is a linear ordinary 
differential equation for ij) vtiQ < y < 1: 


dy'^ 



2 

= R[X + iku{y)] 



Ip — iku” {y)Rip. 


(46) 


The boundary conditions at the free surface are unaffected by a, and after some manipulation involving 
([^ to eliminate the fluid pressure, we obtain three boundary conditions at the free surface: 

-tp'''{l) + {RX+ikR-3k'^)'iP'(l) = 2ikHcot0+——, tp^'{!) =-2H-k^i^{l), iki;{l) = {X+ik)H. (47) 

o 

The no-slip boundary condition on the wall yields 

17(0) = -V''(0) = 0 (48) 


and the responsive flux through the wall becomes the boundary condition 


1/(0) = ik'ip{Q) = —aH. 


(49) 


When fc = 0, we can solve the system (46 )-(H^ for -ip^y) and A analytically, and enumerate the eigenmodes. 


There is a single eigenmode that involves perturbations to the interface height (i.e. H 0), and for this 
eigenmode A = —a at fc = 0. There are also an infinite number of shear eigenmodes which leave the interface 
position unperturbed. These eigenmodes are all stable, and the eigenvalue with the largest real part satisfies 
XR = —{tt/2)‘^, irrespective of a. 

For k ^ 0, we solve the system (46)-(49) numerically. We can formulate the system as a generalised 


eigenvalue problem for ip, H and A, discretise the derivative operators on ip using finite differences or 
Chebyshev polynomials, and solve the resulting generalised eigenvalue problem using standard Matlab 
routines. Alternatively, we can formulate the complete system for the real and imaginary parts of ip and A 
as a boundary value problem in AuTO-07r^ with A as a free parameter, as discussed by Kalliadasis et al.^, 
though this is only useful for tracking a single eigenmode. 

Results for A(fc) are shown in Fig. for the two least stable eigenmodes. In the absence of controls, the 
Navier-Stokes results show a smaller cutoff wavenumber than the Benney and weighted residual results. 
As was the case for the weighted-residual equations, we find that introducing positive a decreases the real 
part of both eigenvalues shown, but has vanishing effect when fc = 0 on all but the least stable eigenmode. 
Furthermore, the critical a computed according to the Benney result (34) is again sufficient to stabilise the 
uniform state against perturbations of all wavelengths. 

In Fig. 1^ we show the critical a required so that 3fi(A) < 0 for all k in the Navier-Stokes equations. This 
is computed in Auto-07p, with the condition that 5ft(A) has both a turning point and a zero at the same 
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FIG. 5. Results for the real part of the two Navier-Stokes eigenvalues A with largest real part, as a function of fc, 
for R = 5, C = 0.05, 6 — 7r/4 and a = 0 (solid lines), and a = as (dashed lines) from (341. 



FIG. 6. Results of an initial value calculation using the weighted residual equations, starting from a non-uniform, 
non-equilibrium state, which evolves without suction until t = 125. For t > 125, we enable feedback controls with 
F = —0.1(h — 1), and the system converges towards the uniform state. 


value of fc. When a < 0.5, this stability boundary is in good agreement with the weighted-residual results, 
with both predicting that the critical Reynolds number is increased substantially, from its uncontrolled value 
of 1.25 to around 50. Beyond this point, the weighted-residual results predict that the critical R should 
continue to increase rapidly with a. However, the Navier-Stokes results show a turning point in R{a), 
followed by a very slow decrease in i? as a is increased. This eventual deviation is not unexpected, given 
the wide range of Reynolds number spanned in this calculation. 


D. Initial value calculations 

Although we have demonstrated that the control parameter a can be chosen to make the uniform state 
linearly stable to perturbations of all wavelengths, it is not necessarily the case that the system will converge 
to the uniform state in nonlinear simulations. In Fig. we show results of an initial value calculation of the 
weighted-residual system, starting from a finite-amplitude state that is neither a steady nor travelling-wave 
solution of the weighted-residual equations. We initially allow this state to evolve without controls, and find 
that the system moves towards a travelling wave state of finite amplitude. We then activate the feedback 
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(a) Benney (b) Weighted-residual 




FIG. 7. Linear stability properties of the uniform state as a fun ction of the control strength a and the displacement 
^ between observer and actuator, with the control scheme (501 for i7 = 5, 9 = njA, C = 0.05. Stability results refer 
to perturbations of all wavelengths. The lowest a is required at a Hnite positive value of The dashed line shows 
the 0{k^) optimiser in the Benney equations: ^ = 2/7/3. 


controls with a = 0.1, which is large enough that the uniform state is linearly stable. After the decay of 
transient behaviour, we observe that the distance of the solution to the desired state decays exponentially 
with respect to time, which is consistent with the expectation that the largest deviation is due to a single 
eigenmode which decays at constant rate. As the imposed injection is proportional to /i — 1, the control 
magnitude also decays exponentially with time, and would be expected to become vanishingly small at late 
times. However, although the amplitude of the applied injection and suction may become very small at late 
times, the feedback control scheme is still required to suppress the growth of small perturbations, and thus 
to ensure the linear stability of the system. 


E. Linear stability for phase-shifted distributed controls 


As an initial step towards designing a more efficient system for feedback control, we can investigate the 
effect of shifting observations relative to actuators, still using a normal mode analysis. We replace the control 
scheme (171 with a scheme based on shifted observers: 


F{x, t) = —a[h{x — t) — 1]. 


(50) 


Here the real parameter ^ is the distance between observer and actuator. Positive ^ means that the observers 
are displaced upstream relative to the position at which the injection is applied. This scheme introduces no 
favoured x locations, and so the eigenmodes can still be written as 

h = 1-\-eheyiY){ikx + \t)-\-0{e^), g = 2/3-b eg exp(fA:a;-|-At). (51) 


We then find 


F =—ae *^^e/iexp(i/ca:-b At). (52) 

We thus simply replace a by aexp(—tfc^) in ( [M| ) and ( [^ to understand the effect of ^ on the eigenvalues of 
the Benney and weighted-residual models respectively. For both models, we can perform a numerical search 
to calculate the boundary of the region in a-^ space where the uniform state is stable to perturbations of 
arbitrary wavelengths, as shown in Fig. for the parameters in this figure, we find that choosing ^ « 2 
has the best stabilising effect in both models, in the sense that a stable uniform state is obtained at the 
lowest value of a. There are some differences between the results for the two models: the effect of positive 
^ is less pronounced in the weighted-residual model than in the Benney calculations, and in fact for the 
weighted-residual model, choosing positive ^ eventually becomes less stabilising as a is increased. 
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In order to understand why the control scheme is most effective when actuators are displaced upstream 
by a finite distance, we can expand the Benney eigenvalue under the assumption that is small, to reach 

3?(A) = -a (^1 + 

To maximise the effect of a, we should choose ^ = 2i?/3, which provides a reasonable estimate of the optimal 
as shown in Fig.[^ This should become a better estimate as R ^ Rq, so that the unstable k move towards 
zero. 


8k^ 

~15 


R- 


5cot6> 


8C 


(53) 


IV. APPLYING CONTROLS VIA POINT ACTUATORS 

A physically important question that we wish to address next is the application of suction controls using 
point actuators, and based on a limited number of observations of the system state. Here we consider only 
behaviour within a spatial period of length L, and only stabilisation of the uniform state. 

We are given the localised actuator functions so that 

M 

F{x,t) = ^ (54) 

m—1 

where the M coefficients bm{t) are to be determined from P discrete observations yp{t) of the interface 
height: 


= [ ^p{x){h{x,t) - l)dx. (55) 

Jo 

We note that the explicit x-dependence of the system that arises from localised actuators and observers 
means that the system is no longer translationally invariant in x, and so linear stability properties of even 
a uniform film in the Navier-Stokes equations cannot be obtained by a normal mode analysis. Instead, we 
derive most of our control strategies using the Benney model, and use the weighted-residual model as a 
black box experiment to represent the additional complexities of the full physical system subject to controls 
derived using a low order model. 

As a starting point, we suppose that the controls are a linear function of the instantaneous observation, 
which is known as a static observation scheme. In the most general form, we can then write 

F = 4'A:$(/i-I). (56) 

Here the operator $ describes observations of the system, 'k represents the shape of the actuators, and 
K is the control operator which we are free to choose based on our knowledge of 4), and the system 
dynamics. We will use M linearly independent actuators and P observations, which are the ranks of 4* and 
4) respectively. In a discretised form, tk and 4> are matrices of size N x M and P x N respectively. The 
matrix K has size M x P, and we may choose all of its entries. Given this form for F, we can compute the 
linear stability of a given steady state by replacing a with — 4/A'4> in the eigenvalue problems described in 

Sec, rrm 


A. Choice of point actuators and observers 

We choose to use M equally-spaced actuators, which are each periodic with period L and locally behave 
as Dirac d-functions, so that 


'km(a:) = <5(x - Xm), Xm = mL/M. (57) 

We similarly use P equally-spaced observer functions, which are displaced upstream by a distance ^ from 
the actuator positions, so that 


4>p(x) = 6{x - Xp), 


Xp = pL/P- C 


(58) 
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FIG. 8. A typical row of the matrix K, or feedback gain, obtained by the LQR algorithm, with 5 equally spaced 
actuators, with shape smoothed according to ( |59[ ) with w = 0.1, and shown by the dotted line here. The cost 
parameter is — 0.1, and for the weighted-residual equation, the same cost weighting is associated with q — 2/3 as 
for h — 1. 


For our numerical calculations, we replace 6{x) in (57) and (58) by the smoothed, periodic function d{x), 
defined by 


d{x) = exp 


cos(27ra;) — 1 


(59) 


One such actuator shape function is plotted in Fig. §for w = 0.1. We normalise the smoothed functions d{x) 
so that each actuator and observer shape function has integral 1 over the interval [0,L], and so d{x) —> S{x) 
as w; —>■ 0. 


B. Proportional control 


If the number of actuators is equal to the number of observers, one of the simplest methods to choose the 
suction/injection prohle is to link each actuator to a neighbouring observer, setting 


bmit) = -aym{t) 


(60) 


where the positive control amplitude a acts in a similar way to the control parameter a in Sec. m In terms 
of the generalised eigenvalue problems, we simply set K = —al. If all actuators, and all observers, are 
equally spaced, the control scheme is specihed entirely by a and the displacement ^ between actuator and 
observer. In Sec. m we considered the continuous analogue of this scheme, with feedback at every point 
proportional to the interface height at that point only. We found that positive a had a stabilising effect on 
the system dynamics according to both long-wave models and also in the Navier-Stokes equations. 

When applying the proportional control scheme (60) with localised observers and actuators, the eigen- 


modes are not sinusoidal, and so we calculate the linear stability properties numerically by solving an 
eigenproblem; note that this calculation only allows for perturbations with wavelength at most L. We Hnd 
that increasing a has a stabilising effect on the uniform Him, and that the value of a required to obtain a 
linearly stable state decreases when increasing the number of actuators M and observers P (see Fig. [^. 

We can also investigate the effect of the displacement ^ between the observers and actuators on the linear 
stability of the uniform state. As discussed in Sec. |III E[ the uniform state is most easily stabilised by 
distributed controls when f « 2, and when R is close to Rq, this best choice for ^ is given by ^ ^ 2i?/3. 
However, the use of localised observers and actuators introduces a natural lengthscale L/M, which is the 
distance between neighbouring observers, and for the analysis in this section, we also have the lengthscale 
L of the imposed periodicity. We can numerically calculate the effect of the displacement ^ on the linear 
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(a) Benney (b) Weighted-residual 




FIG. 9. Stability results as the control amplitude a is varied, with a phase shift ^ between actuator and observer. 
There are M equally-spaced actuators, and P = M equally-spaced observers, each smoothed according to (591 with 


w = 0.1, and results are shown for M — P = 3, 5, 7,9. The largest stable region occurs for M = 9. As is the case for 
distributed actuators (see appendix HIE and Fig. [^, the best stabilisation occurs at a moderate, positive value of 
so that the observers are positioned upstream relative to the actuators. 


stability of the uniform film, subject to the control scheme (601, by solving an eigenvalue problem for each 
^ and a. Fig. shows the stability boundaries in a-^ space for M = P = 3, 5, 7,9. We find that a stable 
state can be obtained at the smallest a when ^ « 2, which is comparable to the results of the calculations 
for distributed controls and observations shown in Fig. despite the additional lengthscales present in the 
system with localised observers and controls. The magnitude of a required to stabilise the uniform state 
generally decreases as M = P is increased, but even for just three actuators, we can stabilise the uniform 
film state by choosing a sufficiently large a with ^ « 3. 


C. Linear-quadratic regulator with full observations 

The control scheme described in the previous subsection only allows each actuator to communicate with 
a single observer. We should be able to obtain better control by allowing data from all observers to be 
combined before determining the actuator amplitudes; we will still consider linear control, but allow all 
entries of the M x P matrix K to be non-zero. This more general scheme can also encompass situations 
where M ^ P. 

The statement that the full system state can be observed is a stringent constraint; for the weighted- 
residual model this requires simultaneous information regarding h{x, t) and q(x^ t), and in the Navier-Stokes 
system, the full system state includes two components of the velocity field along with the interface height. 
Notwithstanding the difficulties of obtaining full observations, if we are somehow able to observe the full 
system state, a variety of algorithms from control theory can be used to compute the controls. Here we 
choose to use the linear-quadratic regulator (LQR) algorithuJ^, which determines K so as to minimise a 
cost functional associated with control amplitudes and the deviation of the system from the flat state. 

The LQR algorithm is designed for the system 

dz 

-— = Jz + 'Fm, u = Kz, (61) 

at 

where z and u are vectors, the matrices J and are given, and we wish to choose the matrix K in order to 
minimise the cost k defined by 


/ {z^Uz + u^Vu) dt, 
Jo 


(62) 


where U and V are given symmetric, positive definite matrices that define the relative cost associated 
with different solution components. A minimiser K of the cost (621 subject to the system (61) is strongly 
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connected to a solution, if it exists, of an algebraic Ricatti equation: 

U + PJ + J^P - P'HV-^^'^P = 0, 


(63) 


definite for all other solutions P, then P is called a minimal solution to 

Furthermore, in Zabczyk^^, it is proved that if the pair (J, 4') is controllable 


(621. 


in which the unknown P is a nonnegative definite matrix. If P is a solution to (63) and P — P is negative 

and K = —V~^'i/^P minimises 

the cost functional 


and U = C'^C, where the pair (J, C) is observable (see Appendix for definitions of controllability and 
observability) then the algebraic Ricatti equation ( [63| has exactly one solution P, and the matrix J — 
is stable. 

For simplicity, we use the following cost functional, in terms of our variables. 


pCO nL 

/ / {^(/i — 1)^ + (1 — ^)P^} dxdt. 

do do 


(64) 


For a given physical system, the control scheme is a function of the single parameter fi € (0,1). The choice 
of AT and the resulting system eigenvalues are dependent on /r, but a stable system should be obtained for 
any 0 < ^ < 1. Row m of the matrix K determines the amplitude of actuator m: 


bm(t) = 


N 


Kmn{hn{t) - 1 ). 


(65) 


Fig. [8] shows one such row, or feedback gain, computed using the LQR algorithm, as implemented using 
the Matlab LQR function, for the Benney and weighted-residual equations. The LQR algorithm gives very 
smooth control input functions for the Benney equation. The largest part of the input function is localised 
slightly upstream of the actuator location when using the full Benney equation (131, or more centrally when 
using the simplihed version ( [T^ . We can insert the Benney controls directly into the weighted-residual 
model, and in fact still obtain a stable state. 

We can also use the LQR algorithm to calculate controls for the weighted-residual model, but the controls 
require observations of both h and q. We also note that the control input functions (Fig. have relatively 
sharp edges near the width of the actuator. The full LQR controls are able to stabilise the uniform state in 
the weighted-residual model, and for our test case the maximum real part of any eigenvalue is —5.62 x 10“^. 
Realistically, we are unlikely to have access to observations of both h and g, and so it would be desirable 
to approximate q from our observations of h using a low order model. The simplest method is to suppose 
that q = 2/3, in effect discarding the control component from q. We find that this yields a linearly stable 
system, but the maximum real part of any eigenvalue is then —5.09 x 10“^, so that convergence towards the 
uniform state would be very slow. We can recover the information regarding the q controls by supposing 
that q = 2h^/3 (which is the leading order term in the long-wave flux ([T^), and so g = 2h. The largest 
growth rate is then —5.64 x 10“^, which is comparable to the growth rate obtained when the flux g can be 
fully observed. 


D. Dynamical observers for a finite number of observations 


For the LQR methodology described above, full observations of the system state are assumed to be 
available. The system is specihed by the interface shape in the Benney equation, but in the weighted- 
residual equations we also require full knowledge of the total down-stream flux at each streamwise location. 
Furthermore, for the Navier-Stokes equation we would need to know the instantaneous velocity at every 
point within the fluid. Such knowledge is unrealistic, and so we now consider the case where the only 
system observations available are those of the interface height, h, at only a finite number of points within 
the periodic domain. In the previous subsection, we showed that if full observations are available, standard 
algorithms, such as LQR, can be used to construct a control matrix K for the instantaneous control scheme 
(56) so that localised actuators can be used to stabilise the uniform state. Alternatively, if distributed 
actuators can be applied, the LQR algorithm can also be used to calculate a control scheme subject to 
localised observers (see discussion in Appendix). However, if there are restrictions on both actuators and 
observers, it is not always possible to construct a control scheme based on ( [5^ so that the uniform state is 
linearly stable. Instead, we turn to a system of dynamical observers, in which both current and historical 
observations are used to determine the controls. 
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The principle of the approach described here is to construct an approximation of the system state which is 
continually corrected based on the observations available. We focus our effort on approximating the coeffi¬ 
cients of those modes which are unstable in the uncontrolled system. We use the dynamic method described 
by Zabczykl^ and applied for the KS equation in Armaou and Christofidesl^, where the predictions evolve 
in time according to our understanding of the linearised system behaviour in the form of its Jacobian matrix 
and the system amplitudes, and the predicted amplitudes are corrected according to our observations. This 
is in contrast to a static observation scheme (56), where the controls are calculated only from the most 
recent set of observations. 

After transformation to Fourier space, we can describe the evolution of a small perturbation h in the 
(simplified) Benney equation (14) by 


dh 

dt 


Jh + F. 


( 66 ) 


In the absence of suction, the system has no preferred positions, and so the eigenvectors of J are Fourier 
modes, and the transformed Jacobian matrix J is diagonal. We reorder the wavenumbers so that the unstable 
eigenmodes of J appear first: 


cm 

dt 


Jh + F = 



(67) 


where the subscripts u and s correspond to unstable and stable modes, respectively. We wish to control to 
the state h = 0. 

To stabilise the zero state of this system, we would ideally leave the stable modes untouched, while 
choosing F to react to the unstable modes. This can be achieved by letting 

F = J/Kh = 

so that 

± (hu\ ^(Ju 0 \ (K\ f^uK 

dt \hj V 0 jJ \hs) \^sK 

The matrix on the right-hand side of the eigenvalue problem is lower triangular by blocks, and the block Jg 
is diagonal. The eigenvalues and eigenvectors corresponding to Jg are thus unchanged by F. 

The remaining task is to stabilise the subsystem 


KK 


( 68 ) 


Ju + JtuK 0 \ I hi 

J>gK Jg \k 


(69) 


dhu 

dt 


— Juhu T ^ud^uh-u 


(70) 


To choose the matrix Ar„, we use the LQR algorithm on the subsystem (70), which has size equal to the 
number of unstable modes, M. However, to apply these controls, we need to approximate z = hu based on 
our observations. We can write our discrete set of observations as y = $(/i — 1), y = ^h = + ^s^s- 


We can obtain a good approximation of ^ by considering a set of ordinary differential equations: 


dz ^ ^ ^ ^ 

= {Ju + J'uKu)z + L{y — y) = yJu + '^uKu ~ 


z + Ly, y = ^iiz. 


(71) 


Here y is the expected set of observations based on our current approximation to the system, and the L{y—y) 
term indicates a correction based on our actual observations. Once we know z, we can set F = '^KuZ. 
However, we still need to choose the matrix L in order that z will converge rapidly to /i„. We define an 
error term: e = hu — z, and after several substitutions we hnd that e is governed by 


de 

dt 


{Ju - L^u)e - L^ghgS = Ye- L^ghg. 


(72) 


To obtain rapid convergence of our estimator z towards the true system state, we need the eigenvalues of 
the matrix Y to have large and negative real part, and we can obtain a suitable matrix L by using the LQR 
algorithm. If these conditions on the eigenvalues of Y are satished, it can be proved that the solution z to 
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FIG. 10. Distance between cnrrent solution and uniform film state as a function of time for M = 5 actuators, with 
P observers; and maximum real part of the eigenvalues of the system (731 as a function of P. The actuator and 
observer shapes are as described by (591, with tn = 0.1 and ^ = 0. The initial condition is h = 1 + 0.3cos(27rx/L) + 
0.1 sin(47ra;/I/), with L — 30. 


equation ( |71[ ) converges exponentially to the true coefficients of h as long as the initial guess is sufficiently 
good. Furthermore, if the real part of these eigenvalues is sufficiently large (in absolute value), then we can 
write Y = Y/e for small e and, by multiplying (721 by e, can obtain a system of equations in the standard 
singularly perturbed fornJ^. This system possesses an exponentially stable fast subsystem (the equation 
for z) and an exponentially stable slow subsystem (the (stabilised) linearised equation for h), which implies 
that the system (73) is exponentially stable. 

We can rewrite the complete system in real space, to determine the behaviour of the nonlinear initial 
value problem: 


ht + qx = 

h? ( 

q{x, t) = — I 2 — 2hx cot 9 + 


C 


SRh^hx 

15 ^ 


F{x,t)=T-^i'K^z, 
dz 
dt 

y = $(/i- 1), 


= J„ + 'ifuKu - L<^>, 


Ly, 


(73a) 

(73b) 

(73c) 

(73d) 

(73e) 


where T is the Fourier transform operator. It can be seen from equations (73c)-(73e) that the feedback 


control F is calculated only from those observations of the system state attainable through the matrix <&. 

It is necessary to alter, and hence approximate, all of the unstable eigenmodes of the system in order to 
stabilise the uniform state, and so the size of 2 must be equal to or greater than the number of unstable 
modes. We expect to achieve better performance as the number of tracked and stabilised modes is increased. 
The number of actuators M need not be equal to the number of observers P, and Fig. shows system 
eigenvalues as P is increased for M = 5 (note that P is odd). We find that choosing P = 7 gives much 
faster convergence than P = 5, but further increases in P have negligible effect on the eigenvalues. However, 
nonlinear initial value simulations of the system (73) benefit from taking P = 9. In Fig. |1I[ we compare 


nonlinear initial value calculations for M = 5, based on P = 5 and on full observations. We find that much 
faster convergence is obtained with full observations. 

The system (73) is presented for the simplified Benney equation (14). The analysis can be extended 


to include cross flow effects present in (13) by left-multiplying 'F by (I — d^Zp) before computing 
The Benney control scheme can be implemented in the weighted-residual equation by simply replacing 
equation (73b) by (16), but we cannot be certain that the resulting system will be linearly stable. For our 


test case, we find that even the linear stability of the uniform state in the weighted-residual equations is 
sensitive to P, with P = 5 stable, but P = 7 unstable. A full analysis of the approximately-controlled 
weighted-residual equation is a topic for future work. 
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FIG. 11. Semi-log plot of the distance between the current and flat states (left) and amplitudes of controls as a 
function of time (right), for M = 5. For the upper row of figures, we use P — 5 observations, while for the lower 
row, we use full knowledge of the interface height h. 


V. CONTROLLING TO NON-UNIFORM SOLUTIONS WITH DISTRIBUTED CONTROLS 


Feedback controls of the form F = —a(h — H) can also be used to drive the system towards non-uniform 
states, by setting the target state H to be spatially varying. We would like to know whether the state h = H 
is always reached, and whether this state is stable. As F = 0 when h = F[, the system can only remain in 
this state if = iJ is an exact solution of the equations in the absence of suction. Small perturbations about 
F[ are always affected by the feedback controls, and so a will change the linear stability properties of the 
state F[. We will discuss the system dynamics when FI is not an exact solution of the governing equations 
in Sec. |V C[ The extension of the localised actuator control scheme developed in the previous section to 
stabilise a non-uniform state is a non-trivial task, as discussed in Sec. VD[ and is left for future work. 


A. Travelling waves 


The long-wave systems support non-uniform travelling wave solutions, of the form h = H{x — Ut), where 
U is the propagation speed. Travelling waves undergo bifurcations (see e.g. Oron and Gottlieb!^, and 
may be stable or unstable in the corresponding moving frame. It is important to note that the shapes 
and bifurcation structure of travelling waves differ between the models. If the target state H is an exact 
travelling wave solution to the equations in the absence of suction, then the state H is also a travelling 
wave solution to the same equations with F = —a\h — ^^(C)], and thus the application of controls affects 
the stability but not the shape or speed of the targeted travelling wave. 

Fig. 12 a) shows an unstable travelling wave solution to the Benney equation. For simplicity, we limit 
perturbations to those periodic with the same spatial period as the travelling wave. In order to compute 
the stability of travelling waves, we transform to the frame moving at speed U, and then identify x with ((. 
For the Benney equation, the generalised eigenvalue problem (24) becomes 


Xh = {Udi: - dc^Zh -[I- dcZF]a} h. 


(74) 


We note that if Zp = 0, which is the case for the simplified Benney equation (141, then the eigenvalues A 
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FIG. 12. (a) A travelling wave solution to the Benney equation, for i? = 2, 0 = 7r/4, (7 = 0.05, U = 2.82. (b) The 
real part of the seven complex eigenvalues with largest real part, as a is increased. Real eigenvalues are shown by 
red dashed lines, while black solid lines indicate the real part of complex conjugate pairs. Neutral stability occurs 
at a = 0.0434. (c) Results from nonlinear initial value calculations, starting from close to a uniform film, controlling 
towards the solution shown in (a), for a = 0, a = 0.05, a = 0.1, a = 0.15. Convergence to H is only achieved in the 
two latter cases. 


are shifted by —a, and the eigenvectors of the system are unchanged from those in the absence of controls. 
However, if we are using the standard Benney equation (131 or the weighted-residual system, the effect of a 
on the eigenvalues is more complicated than a simple shift, and we solve the eigenvalue problem numerically 
to determine the effect of a on the linear stability properties of the non-uniform travelling waves. 

Fig. |12[ b) shows the real part of the seven most unstable eigenmodes as a function of a when considering 
the linear stability of the travelling wave shown in Fig. l^a). When 0 = 0, this travelling wave is unstable, 
with one eigenmode with a positive real eigenvalue. There are two eigenmodes with eigenvalue zero; one 
corresponds to varying the mean film thickness, and the other to translational displacement of the travelling 
wave. The real part of the most unstable eigenvalue decreases with a, until it collides with another eigen¬ 
value while still in the right half plane. These two eigenvalues then form a complex conjugate pair, which 
eventually crosses the imaginary axis with a finite imaginary part, stabilising the system for a > 0.0434. 
This stabilisation occurs via a Hopf bifurcation, and so we would expect to observe small-amplitude limit 
cycles for a just below the critical value. However, linear stability alone does not mean that the travelling 
wave is necessarily an attractor when starting from the uniform state, and indeed initial value calculations 
starting from the uniform film state, as plotted in Fig. 12 c), do not reach the desired travelling wave when 
a — 0.05. However, the system successfully converges to the desired travelling wave when a = 0.1, and 
converges more rapidly when a = 0.15. We note that even when the system is converging to the travelling 
wave, the solution norms show evidence of decaying oscillations. 


B. Non-uniform steady states 

We do not know of any non-uniform steady states to the Benney, weighted-residual or Navier-Stokes 
equations for flow down a planar, unpatterned wall in the absence of suction; instead structures are swept 
downstream by the underlying flow. However, as discussed in Thompson, Tseluiko, and Papageorgiou^, 
the application of steady non-zero suction gives rise to non-uniform steady states, with their own bifurcation 
structure and stability properties. Moreover, we can often choose the applied steady suction in order to 
make a given interface shape into a steady solution of the equations. 

In order to include both steady suction and feedback, we use an extension of the controls: 

F =—a[h{x,t) — H{x)] + S{x). (75) 

Here a is the control parameter, and S{x') is the steady component of F that we are free to specify. If a = 0, 
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FIG. 13. Steady flux q and suction S for the steady state (761. The solid and dashed lines correspond to Benney 
and weighted-residual results, respectively. 


S{x) must have zero mean to prevent growth in fluid mass, and thus to allow steady solutions. 
We choose the following non-uniform steady state as the target state for our calculations: 


27r 


H{x) = 1 -I- 0.3 cos — I -I- 0.2 sin — ) -I- 0.2 sin — ) , L = 30 




Gtt 


(76) 


shown in Fig. 13 and set R = 5, C = 0.05, 6 = tt/A. For these parameters, the uniform film state is 
unstable. The state h = H is not a steady solution of the equations when ^ = 0, but we can calculate S{x) 
to make it so. 

For h = H tohe a. steady solution of the Benney equation, we have 


F = S = Qx 


and 


g = — f 2 - 2Hx cot 0 + 


R 




We can rearrange these two equations to obtain a single equation for S = F: 


S- 


/2RH^S 
V 3 




2 — 2Hx cot 9 + 


Hx 


C 


SRH^Hx 

15 


(77) 


(78) 


(79) 


subject to periodic boundary conditions on S{x). The right hand side of (79) is known, and the left hand 
side is linear in S{x). There is therefore a unique solution for S{x), given H{x), in the Benney model, and 
the equation has a solution for each smooth, non-zero H. We note that the task of finding a suction profile 
to enable a particular steady solution is related to inverse topography problems, in which the bottom profile 
is computed from observations of the interface heightP or surface velocitj^^. 

Perhaps unsurprisingly, the linearity with respect to S obtained in ( [7^ does not apply in the weighted- 
residual model. Instead we must solve 


F = S = qx 


and 


q = —— ( 2 — 2F[x cot 9 - 


Hx 


C 


R 


flSq^Hx MHqqx , HqF\ 


\ 35 


35 




(80) 


(81) 


again subject to periodic boundary conditions on S and q. We can use F = qx to rewrite the second equation 
as an equation for q alone: 




q = 


2 — 2Hx cot 9 + 


H. 


XXX 


R 


(ISq'^Hx 

V 35 


27 Hqqx 
35 


(82) 
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FIG. 14. Illustration of a transcritical bifurcation that occurs when controlling to an unstable steady state (with 
minimum value 0.1375), that has just one positive eigenvalue. Exchange of stability occurs through a transcritical 
bifurcation at a = 1.92, necessitating the existence of another solution branch, which here connects to a stable 
steady solution for the same S' at a = 0. The second solution branch only persists slightly beyond the transcritical 
bifurcation, here diverging through the minimum layer height vanishing at a finite value of a. The parameters here 
are R = G, C = 0.05, 9 = n/A, S = 0.7 cos(27rx/10), which matches the bifurcation structure for a = 0 shown in Fig. 
3 of Thompson, Tseluiko, and Papageorgiou^^S. 


but this is nonlinear in the unknown 5 , and so we cannot guarantee existence or uniqueness of solutions. 
However, solutions should still exist for H near to 1, and for the non-uniform state (76 1 , we obtain the 
solution shown in Fig. |13| 

With the appropriate S for the corresponding model, as shown in Fig. |13[ numerical solutions of the 
discretised eigenvalue problems described in Sec. HE show that the steady state (76 1 is stable for a > 1.32 
in the Benney model, and a > 1.39 in the weighted residual model. In both cases, the exchange of stability 
occurs via a Hopf bifurcation, so that below the critical value of a, we would again expect to observe 
time-periodic limit cycles. 

A second mechanism for exchange of stability involves real eigenvalues passing through zero. In Fig. |14[ 
we choose a steady flux S{x) which is knowrPS to correspond to two solutions H{x) when a = 0, one of 
which {Hi) is stable, the other {H 2 ) is unstable with one positive real eigenvalue, and show the results of 
controlling towards the latter, unstable state, H 2 . Each steady state at a = 0 gives rise to a solution branch 
for a > 0. The target state H 2 is always a solution, and is stable for a > 1.92. For a < 1.92, H 2 has a 
single eigenvalue with positive real part, and this eigenvalue is exactly zero at a = 1.92. The exchange of 
stability via a real eigenvalue passing through zero corresponds to a transcritical bifurcation, and implies 
the local existence of a second solution branch, which here is the branch that connects back to Hi dX a = 0 . 
The second branch diverges as a increases beyond 1.92, here by the minimum film height tending to zero. 


C. Controlling towards non-solutions 


In the previous subsection, we assumed that the target state H is an exact solution of the equations, 
so that the system will remain at h = 77 if it ever reaches it, and the main questions surround linear 
stability, which can be directly modified by linear feedback controls. However, in practice, the target state is 
highly unlikely to be an exact solution, due to discretisation error, imperfectly known parameters, or, more 
interestingly for our purposes, discrepancies which arise due to calculating travelling waves or the steady flux 
S according to a low-order model which only approximates the true system. We now investigate robustness 
to model choice by analysing the behaviour of the system when feedback controls are applied towards a state 
which is not a solution to the governing equations, and so can never be more than transiently achieved. 

We suppose that the system reaches an equilibrium state 77*, which will depend on the target state 77, 
the feedback control strength a, any patterning imposed on the system via S, and the parameters of the 
uncontrolled system. We usually have a nonlinear system to solve for 77*, which need not have unique 






solutions. In the Benney equation, the steady state H* must satisfy 

F = -a[H^ -H] + S, F = q^ 
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(83) 


and 


]Lr*3 / Tj* 

q = — h-2H:cote + ^ 


R 


8H*^m 2H*'^F 
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(84) 


These equations are nonlinear in H*, and can have zero, one, or more solutions. Fig. |15| shows steady 


solutions to the nonlinear Benney system (83l-(84), and also the corresponding weighted-residual system, 
for the case S = 0, with H given by the large-amplitude, non-uniform state (76). We find that for both 
models, the numerical solutions for H* tend towards H as a increases, and our linear stability calculations 
confirm that the states H* are stable at large a. However, the value of a at which steady states become 
stable, and also the extent to which the steady states deviate from H at a given a, are dependent on the 
choice of model. 

The linearity of the control scheme means that suction can be interpreted as feedback controls towards 
the equilibrium state H*, which is itself dependent on a and the original target state H. We define S* = 
—a{H* — H) + S, so that for general h, we can write 


F = -a{h -H) + S = -a{h - H*) + S* 


(85) 


As a result, the system is indistinguishable from controlling to the state H* with feedback control parameter 
a and steady suction component S*. 

For large a, we can find a simple asymptotic solution for the steady state JI*. If the system tends towards 
a bounded steady state as a increases, then F must remain bounded, and so the interface shape H* must 
tend towards H. Also, F tends towards Sq{x), which is defined to be the steady flux required to make the 
desired state FI a steady solution of the equations. Thus, without regard to the model details, but assuming 
only that a bounded steady state F[* exists for large a, we find that this state behaves as 


H* = H + 


S-Sa 


a 


o\ 4 


( 86 ) 


The function Sq, and subsequent terms in the expansion, will depend on the details of the model, but in 
general we can move the equilibrium state H* closer to the desired state F[ by increasing a. In Fig. we 
show a system where there are two steady states for the same parameters, and controls are applied towards 
one of these states. However, one of the solution branches disappears at a finite a, so that for sufficiently 
large a, the only steady state remaining is the one described by (86). More generally, branch divergence 


means that unwanted solution branches can be eliminated by increasing the control amplitude. 


D. Control towards non-uniform states with point actuators 

In Sec. m we considered control schemes based on localised observers and actuators that remain fixed 
in the laboratory frame, and showed that these schemes can be used to stabilise the uniform film state. 
We then showed in Sec. |^that distributed control schemes can be used to stabilise non-uniform travelling 
waves, and to create and stabilise non-uniform steady states. However, the extension of the point-actuator 
control schemes to non-uniform travelling waves and non-uniform steady states faces significant difficulties. 

Travelling waves are steady with respect to a moving coordinate (^ = x — Ut, and can be written as 
h = H{(). However, if the observers and actuators are fixed in the laboratory frame, then these move 
relative to the travelling wave to be controlled. To calculate linear stability, we first transform to the 
moving frame, so that the base state h = H{Q is a steady solution of the controlled equations. However, 
the evolution of small perturbations is subject to the spatial structure of the control scheme, which in this 
moving frame is also time-dependent. If the control scheme is spatially periodic in the laboratory frame, 
then it is both spatially and temporally periodic in the travelling frame, and we must use Floquet multipliers 
with respect to time to obtain eigenvectors. For a general control scheme, this requires the computation of 
eigenfunctions that are explicitly dependent on both space and time within periodic boundary conditions, 
which is beyond the scope of the present study. We note however that Gomes, Papageorgiou, and Pavliotis® 
showed that for the KS equation, localised controls derived for a uniform state could be used to stabilise 
travelling wave solutions in cases where the wave frame moves relative to the controls. 
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FIG. 15. Steady solutions to the Benney equation and weighted residual equations, controlled towards h = H (shown 
in bold) using the control scheme (751, for a = 0, 0.125,0.25, 0.5,1, 2, and S' = 0. Here 7? = 5, C = 0.05 and 6 = 7r/4. 
Dashed lines indicate unstable solutions. 


Non-uniform steady interface shapes H(x) require a non-zero suction profile So{x) in order to be steady 
solutions of the governing equations. However, if suction must be delivered through a linear combination of 
M localised actuator shapes, it is very unlikely that the exact profile So{x) can be achieved. Thus we will 
no longer obtain the result that h —>■ H(x) when strong controls are applied. It is easy to imagine situations 
where the interface shape appears to be close to the desired state when viewed through localised observers, 
while diverging significantly at other positions, and so we leave the analysis of this system to future work. 


VI. CONCLUSION 

In this paper, we analysed the effect of feedback control on the dynamics of a thin film flowing down 
an inclined plane. Feedback was applied through injection and suction through the planar wall, with the 
required injection/suction profile determined in response to observations of the position of the air-fluid 
interface. We used long wave models based on the Benney and weighted-residual methodologies to describe 
the effect of suction and injection on the system dynamics. We note that suction is the only mechanism by 
which the net system mass can be modified, and so suction controls are the only way in which perturbations 
of infinite wavelength can be made better than neutrally stable. 

The simplest control scheme is to suppose that the suction profile is locally proportional to the deviation 
of the interface profile from the desired state, so that fluid is injected where the film is particularly thin, 
and removed from thicker regions. We used a linear stability analysis to show that this simple control 
scheme, governed only by the constant of proportionality a, has a stabilising effect on the uniform film state 
for positive a in both Benney and weighted-residual models, and also in the Navier-Stokes equations. We 
calculated the critical value of a needed to stabilise the uniform state to perturbations of all wavelengths, and 
showed that the control scheme can significantly increase the Reynolds number for the onset of instability. 

The analysis summarised so far is for distributed controls, but we also studied a more realistic scenario 
by supposing that injection/suction can be delivered only via a small number of localised actuators, corre¬ 
sponding for example to slots in the planar wall. Likewise, we should base our control scheme on a limited 
number of observations of the system state. The control system requiring the least amount of communica¬ 
tion between actuators and observers occurs when each actuator is connected to only one observer, and the 
applied suction is proportional to the deviation of the observation from the desired value. For equally spaced 
actuators, our numerical calculations show that this singly-connected control scheme has a stabilising effect 
on the uniform film state in both long-wave models. The uniform state becomes more stable as the number 
of observers and actuators is increased. We investigated the effect of displacing the observer relative to its 
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linked actuator, and found that the observer should ideally be positioned slightly upstream of the actuator 
to obtain the best stabilisation. Displacement between observers and actuators can also be incorporated in 
the fully distributed case, and we again Hnd that the most efficient stabilisation occurs when the observer 
is slightly upstream of the actuator. 

In principle, we should be able to obtain better system performance by using all available observations to 
compute the feedback controls. If the entire system state is observable, we can use standard algorithms from 
control theory to decide the control inputs according to various objectives. For example, we used the linear 
quadratic regulator (LQR) algorithm to minimise a cost functional defined in terms of the deviation of the 
film from the flat state and the actuator amplitudes from zero. The use of point actuators means the system 
is not translationally invariant, and so the linear stability of the Navier-Stokes equations can no longer be 
studied by a normal mode analysis. However, we found that controls calculated using the LQR algorithm 
for the Benney equation were able to stabilise the uniform state in both the Benney and weighted-residual 
systems. 

For the case where only a small number of observations are available, controls developed under the 
assumption of full observations can still be implemented by using dynamical observers, and we exploited 
this strategy to control the Benney system. In this scheme, the Benney system is augmented by a system 
of ordinary differential equations to create an evolving approximation of the magnitudes of the unstable 
eigenmodes, which evolves according to our understanding of the underlying system, with corrections due 
to the available observations. Our stability and initial value calculations confirm that this approach does 
indeed stabilise the uniform state in the Benney system. For our test case, we found that increasing the 
number of observations above the number of unstable modes initially yields a significant increase in the 
overall convergence rate, but further increases have negligible effect. 

To test the robustness of the dynamical observer scheme in a proxy physical setting, we inserted the 
Benney control scheme into the weighted-residual equation. We found that the uniform state was sometimes 
stable, but this depended sensitively and non-monotonically on the number of observations used to calculate 
the controls. The eigenvalues of the Benney and weighted residual equations behave differently, and so we 
might expect that the dynamic approximations converge poorly to the true state. However, at least for 
stabilising the uniform state, we have the option of using the singly-connected control scheme with discrete 
actuators and observers, which behaves similarly in both long-wave models, and so depends relatively weakly 
on model details. 

The thin-film systems can support non-uniform travelling waves, which propagate down the slope at 
constant speed. These may be stable or unstable; and we found that the locally proportional controls can 
be used to stabilise unstable travelling waves. The total magnitude of the imposed suction will vanish as 
the target state is approached if it is an exact solution of the equations, so controls can be in principle 
be used to physically verify the shape of unstable states. If a steady suction prohle is applied, the system 
can support non-uniform steady state^^. These steady states have their own bifurcation structure, can be 
stable or unstable, and have a more complicated internal flow than that for a film of uniform thickness. If 
the suction prohle corresponding to a desired steady interface shape is known exactly, we showed that the 
feedback control scheme can be used to stabilise the steady state in a similar manner to that for stabilising 
travelling waves. 

The shape and speed of travelling waves, and the suction prohle corresponding to steady states, differs 
between the two long-wave models here, and likely also the Navier-Stokes equations. It is therefore un¬ 
reasonable to assume that the target state is an exact solution of the equations. However, we hnd that if 
controls are applied with large positive a towards an arbitrary state, the system will both move towards 
that state and become stable as a is increased, irrespective of the model used. 

In some ways it is unsurprising that simple control schemes can be used to linearly stabilise the uniform 
state in the Benney equation, as the linear operator is similar to th at for the Kuramoto-Sivashinsky (KS) 
equation, where control schemes have been rigorously derivecP^^. However, the KS results provide no 
guarantee on the nonlinear behaviour, or on system dynamics away from long-wave limits, and so our 
nonlinear initial value calculations and linear stability calculations in the Navier-Stokes equations provide 
meaningful tests on the use of feedback control. It would be interesting to investigate the effect of linear 
suction controls on nonlinear stability and blow up phenomena in the Benney equation. 

The motivation for this work was to act as a first step towards experimental implementation of feedback 
controls in thin film flow, and indeed the results are very promising, but there are a number of physically- 
motivated questions still to be addressed. The analysis in this paper concerns only 2-D flows, but in 3-D 
simulations or experiments there are more choices regarding the configurations of actuators and observers 
which may be arranged as points or slots, and we would also need to consider edge effects. In Sec. |HI[ 
we analysed linear stability of a uniform film within an infinite domain, while we later applied periodic 
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boundary conditions for our analysis of control strategies using localised actuators and observers, and for 
control towards non-uniform states. However, experiments are actually performed on a wall of finite length, 
with an inlet (at which perturbations can be applied) and outlet. Future work could include exploring the 
control strategies described in this paper with more realistic boundary conditions, assessing the effect of 
noise, and also incorporation of restrictions on the control scheme to reflect latency in flow visualization, 
data processing, and the actual physical setup by which suction and injection is applied. 
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Appendix A: Controllability and detectability 


Here we state some basic definitions from control theory; further details can be found in 
consider the linear system 


Zabczykl^. We 


= Az + Bu, y = Cz, 


(Al) 


where A, B and C are N x N, N x M and M x P matrices, respectively. We will say that a matrix A is 
stable if all its eigenvalues have negative real part. 

We will call the system (Al), or the pair (A, B), controllable if there exists a matrix K such that A + BK 
is stable. If the system is controllable, we can always obtain the state z* by taking u = K{z — z*), regardless 
of initial conditions. Similarly, we say that system (Al), or the pair (A, C), is detectable if there exists a 
matrix L such that A -|- LC is stable. If the pair (A, C) is detectable, then (A^, C^) is controllable. 

The Kalman Rank condition gives a necessary and sufficient condition on A and B for controllability, and 
therefore detectability. This condition states that the system ( |A1[ ) is controllable if and only if rank[A|il] = 
N, where 


[A\B] = [B AB A^B ■ ■ ■ A^-^B] 

is a N X N'^ matrix obtained by writing consecutively the columns of the matrices A^~^B, n = 1,..., N. 

The natural choice when constructing controls based on the observations y would be to choose a matrix 
K such that the matrix A -|- BKC is stable. Controls that can be written in the form u = Ky are called 
static output feedback controls. However, for nontrivial B and C, it is not possible, in general, to construct a 
matrix K so that A -|- BKC is stable. This difficulty motivates the construction of the dynamical observers 
presented in Sec. |IVD| 
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